* FigureA2A3.do  06/01/23
* Sensitivity analysis for Figure 2
* Reported in the on-line appendix as Figures A2 and A3

do figure2

gen rxshalt = deathsrx/deathsall - polysh
label var rxshalt "Rx-only deaths as share of all opioid deaths"
gen xregpool=.
gen xregpoolse=.
replace xregpooln=.
gen lnxrevpool=.
gen lnxrevpoolse=.
gen lnxrevpooln=.
gen lnxaltpool=.
gen lnxaltpoolse=.
gen lnxaltpooln=.
sum deathrtall [aw=pop]
scalar avgdeathrt = r(mean)
foreach yr of numlist 1999/2021 {
  reg deathrtall rxsh i.female##i.agenum if year==`yr' [aw=cumdeaths]      // level specification
  replace xregpool = _b[rxsh]/avgdeathrt if year==`yr'
  replace xregpoolse = _se[rxsh]/avgdeathrt if year==`yr'
  replace xregpooln = e(N) if year==`yr'
  reg rxsh lndeathrtall i.female##i.agenum if year==`yr' [aw=cumdeaths]    // reverse regression
  replace lnxrevpool = _b[lndeathrtall] if year==`yr'
  replace lnxrevpoolse = _se[lndeathrtall] if year==`yr'
  replace lnxrevpooln = e(N) if year==`yr'
  reg lndeathrtall rxshalt i.female##i.agenum if year==`yr' [aw=cumdeaths] // leave mixed consumption out of Rx
  replace lnxaltpool = _b[rxshalt] if year==`yr'
  replace lnxaltpoolse = _se[rxshalt] if year==`yr'
  replace lnxaltpooln = e(N) if year==`yr'
}

gen xregu = xregpool + (1.96*xregpoolse)
gen xregl = xregpool - (1.96*xregpoolse)
gen xrevu = lnxrevpool + (1.96*lnxrevpoolse)
gen xrevl = lnxrevpool - (1.96*lnxrevpoolse)
gen xaltu = lnxaltpool + (1.96*lnxaltpoolse)
gen xaltl = lnxaltpool - (1.96*lnxaltpoolse)


#delimit ;
twoway (rcap xregu xregl year if female==1 & agenum==1, lcolor(black))
  (scatter xregpool year if female==1 & agenum==1, msymbol(S) mcolor(black))
  (line xregpool year if female==1 & agenum==1, lcolor(black))
  (rcap xaltu xaltl year if female==1 & agenum==1, lcolor(blue))
  (scatter lnxaltpool year if female==1 & agenum==1, msymbol(S) mcolor(blue))
  (line lnxaltpool year if female==1 & agenum==1, lcolor(blue))
  (line lnxregpool year if female==1 & agenum==1, lcolor(gray))
  (line zeros year if female==1 & agenum==1, lcolor(gray) lpattern(dash)),
  title("Figure A2.  The cross-area relation between opioid death" "rates and their Rx-Im composition, by year") subtitle("Alternative specifications")
  legend(order(2 "Point est: level spec" 1 "95% confidence interval"
  5 "Point est: alt share" 4 "95% confidence interval"
  7 "Point est: from Fig 2") region(lcolor(white)) size(small))
  graphregion(color(white)) ylabel(,nogrid) xtitle("") ytitle("Cross-area regression coefficient")
  note("Notes: Observations are age group by gender by Census division, weighted by cumulative opioid deaths 1999-2021."
  "Figure 2's series (log mortality rate spec) is shown together with results for a level spec and an alternative share spec." "The horizontal axis is shown as a dashed line.  Figure 2's confidence interval is not shown."
  "The level spec uses the mortality rate level rather than log mortality rate, with coefficient rescaled by mean mortality for comparability."
  "The alternative Rx share is the ratio of Rx deaths, excluding any deaths coded as both Rx and Im, to all opioid deaths."
  "Source: CDC Wonder.", size(vsmall))
  name(figurea2, replace);
#delimit cr

#delimit ;
twoway (rcap xrevu xrevl year if female==1 & agenum==1, lcolor(gray))
  (scatter lnxrevpool year if female==1 & agenum==1, msymbol(S) mcolor(black))
  (line lnxrevpool year if female==1 & agenum==1, lcolor(black))
  (line zeros year if female==1 & agenum==1, lcolor(gray) lpattern(dash)),
  title("Figure A3.  The cross-area relation between opioid death" "rates and their Rx-Im composition, by year") subtitle("Reverse regression specification")
  legend(order(2 "Point estimate" 1 "95% confidence interval") region(lcolor(white)) size(small))
  graphregion(color(white)) ylabel(,nogrid) xtitle("") ytitle("Cross-area regression coefficient")
  note("Notes: Observations are age group by gender by Census division, weighted by cumulative opioid deaths 1999-2021."
  "Each year, the Rx share of opioid deaths by is regressed on gender-age group interactions and log opioid deaths per 100K."
  "The horizontal axis is shown as a dashed line." "Source: CDC Wonder.", size(vsmall));
#delimit cr

